{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Auxiliary Mesh Data Structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We discuss ways to extract the combinatorial structure of a triangulation by using `elem` array only. Auxiliary data structure includes: \n",
    "- 2D: `edge, elem2edge, edge2elem, neighbor, bdEdge`\n",
    "- 3D: `face, elem2face, face2elem, neighbor, bdFace`\n",
    "\n",
    "They are wrapped into a mesh structure `T` and generated by\n",
    "```\n",
    "    T = auxstructure(elem);  % 2-D\n",
    "    T = auxstructure3(elem); % 3-D\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The auxiliary data structure can be constructed by *sparse matrixlization* efficiently; see [Auxiliary Mesh Data Structure](auxstructuredoc.pdf) for detailed explanation. In the following, we present two examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Edge\n",
    "We first complete the 2-D simplicial complex represented by `elem` by constructing 1-dimensional simplices, i.e., edges of the triangulation. We use `edge(1:NE,1:2)` to store indices of the starting and ending points of edges. The column is sorted in a way such that for the k-th edge, `edge(k,1)<edge(k,2)`. The following code will generate an `edge` matrix. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gkHFRQu7wrx4AAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1TZXAtMjAxOCAxNDoyMDo0NufgpwsAACAA\nSURBVHic7d1/eFx1nejxz7TlRylNDwxUEhh6+NG0GNChAbyJYib7+Ah6G+xdMeB1NZNFi+h6lxph\n/YUJ+3h9HpcNqOteoXvtTr2uhcCFaHO1spacgEwoNiWUltJA21OmJPwamAb6u83cPw4MQ5Im8+Mk\nc873+349/hFmvpnnC37n+845Z34E0um0AABQajNKPQEAAEQIEgDAIwgSAMATCBIAwBMIEgDAEwgS\nAMATCBIAwBMIEgDAEwgSAMATCBIAwBMIEgDAEwgSAMATCBIAwBMIEgDAEwgSAMATZpV6AlCWbduW\nZe3evduyLNM06+rqotFoqSeF6cMCQL4CfEEf3GXb9urVqy3LsixLRIJGeaVZnUwNDdh9ImKaZiQS\naWpqikQipZ0npggLAAUjSHCNZVnNzc22bQeN8ppwg4g0RJZn7k2mBgfsTdvtjc7eZJpma2srfzKr\nhAWAIhEkuMOyrPr6+qWR5bXhpUGjYuLBydRgvL9ri/1wNBptbW2dnhliSrEAULyZbW1tpZ4DfM/Z\njFqi99SGG045ee6k4085ee4is3qWnNK5LrZ3717O3vgdCwCuIEgoVmYzqjSr8/rF0FmVbEkKYAHA\nLQQJRSl4M3KwJfkdCwAuIkgoXJGbkSOzJRmGEQ6HXZwephoLAO7iRQ0oXCAQiC5rqwkvLf6hevu7\nrP7Yrl27in8oTBsWANzFG2NRoFgsJiKV5pJJR255/vH/9+gvx95+3dUt5tlVzs+V5pJYZ5tlWZy3\n8YvcF8BDf/r5wO5N2f93Z7z+5uCqh25Lp9PXfvLvnTfSsgB0RpBQoJ6eHhGZ9AW+IjK8742dic1j\nbz9w8O3Mz0GjotKs7unpYT/yixwXgP3S1vUb1hw5cij7/27HyMjIqge/vyOxWUROnW2wAECQUCDL\nsnI8V3PhueGmZe+91+SV13ev+3Nszinzyuefnz0saJTHYjHeleIXky6A9U/85vnd/c8MPHb02JFx\nB/zhsVU7sv5SYQGAD1dFgWzbXmRelsvI+aeHasMNmf8NvrojEJhxw1//0Jh7Zvaw2nCDbdtTMldM\ngUkXQLy/66ltjxyvRrv2bOnq+beLF340cwsLAAQJhcj9+sEo23Zs2Dzw2JXVy6ourBl1V9AoFxHn\nA9DgcbksgM9dteKr193x1evuGHvXocP7Vz1425zZZV+65rZ3bgoEWADglN2UsyzLOduuEsuyKs3q\nXC4gZRsZGen4452BQOATNX8z9l7nMtLtt9+u3n8u9eSyABafd/nx7rpvXfurbyRuuv6fy049PXOj\nwgugrq6Oa2O5IEhTrr6+PtIy+mjA7/p3bL0o+PF8f+vPmzoHX93xocorPxA8d9wBQaO8f8ejUn2o\n6Al6RSox3N+x1awJmbXnlHoubipsATie2vbI45t+W3vpNeHFkXR6JPsu9RaAiFjtvSLCG2xyQZCm\nQ6SlttRTcFkqMTzwcF++vxXv/52IhBdHjjdgwO4zPxlS5j9XKjEcu/Y+p0bK/Es5ClsADuvJ+0Xk\nzb0v333frSLvbNNrfv/jSxZ+TLEFICKdN68LN1b1d2wt9UT8gWtIKES4sSqZGsrrV44cOfTi0HMi\nsuj4Z3KSqSGzNlTs5LzBqdGyu65W7NjIUcACyEhLWkS27XzyqW2PPLWt27lx244Nuwe3qbQARKTz\n5nUisuwnV5d6Ir7BERIKYYTmiciA3Zf7Z8bYg1uPHTtqlM0/47TxLzz09neJiFmjwn6UVaOQ3Zso\n9XTcV8ACyGiILI9c3uj8nE6PrLz/2yJy/adueXt/SlRZAEKNCkKQUAgjVGbWhOL9a3Pfj1548Wl5\n96V049pubww3VhmhMnemWDrZNSr1XKZKAQsgY+GC916bl7mG9MELa/7w2Co1FoBQo0Jxyg4FMmvP\ncb6UOkcvvNgvIqfNnX+8AXk9mmfpUCNHvgtgUmosAKFGReAICQUya0JWe2/uJ22+8YWfTjwgmRpq\naLzSjamVjD41knwWwD1tG493VyAwI3OvAgtAqFFxOEJCgYzQPCNUFu9f68qjvXMByc/7uFY1EhbA\neKhRkQgSCmSEyiIttQN231prZZEPNWD3xTrbog80ujKxktCtRsICGIMaFY9TdihcuLFKRB5rWysi\nDZHlhT3IgN3XHrsx+kCjf7dyDWvkYAFkUCNXECQUpcgtSYHNSNsaOVgAQo3cQ5BQrIK3JAU2I81r\n5NB5AQg1chVBgguyt6RFZvWkL7tKpgbj/V1d1kpfb0bUKEPPBSDUyG0ECe4IN1YZ55R1rri3y1oZ\nNMorzepF5mXZX+CWTA0O2JteTw12vXsN/OYNX/HvuyCp0Si6LQChRlMgwGfQTrVAINA22FLqWUyT\nVGI4ldhr9ybs+B7nI3OcvSmZGnLe9miEysKNVWZNyNf7eF41strjouIH7I5LkwUgedaoraKdnTYX\nHCHBTUaozAiVmbUhaZFUYtjuTdjxxAHZVREqq5AaNTZljo0moMMCEI6NpgxBwlQxQmXhUJVzdUEZ\n1Ch3Si4AoUZTiTfGArmiRqBGU4ogATmhRqBGU40gAZOjRqBG04AgAZOgRqBG04MgAROhRqBG04Yg\nAcdFjUCNphNBAsZHjUCNphlBAsZBjUCNph9BAkajRqBGJUGQgPehRqBGpUKQgPdQI1CjEiJIwDuo\nEahRaREkQIQagRp5AEECqBGokScQJOiOGoEaeQRBgtaoEaiRdxAk6IsagRp5CkGCpqgRqJHXECTo\niBqBGnkQQYJ2qBGokTcRJOiFGoEaeRZBgkaoEaiRlxEk6IIagRp5HEGCFqgRqJH3ESSojxqBGvkC\nQYLiqBGokV8QJKiMGoEa+QhBgrKoEaiRvxAkqIkagRr5DkGCgqgRqJEfESSohhqBGvkUQYJSqBGo\nkX8RJKiDGoEa+RpBgiKoEaiR3xEkqIAagRopgCDB96gRqJEaCBL8jRqBGimDIMHHqBGokUoIEvyK\nGoEaKYYgwZeoEaiReggS/IcagRopiSDBZ6gRqJGqCBL8hBqBGimMIME3qBGokdoIEvyBGoEaKY8g\nwQeoEaiRDggSvI4agRppgiDB06gRqJE+CBK8ixqBGmmFIMGjqBGokW4IEryIGoEaaYggwXOoEaiR\nnggSvIUagRppa1apJ6As27Yty9q9e7eIdN68zqwNhRurSj0pr6NGUKlGqcSw3ZtIJfaKSHNzc11d\nXTQaLfWkPC2QTqdLPQel2La9evVqy7IsyxKRoFFeaVYnU0MDdp+IGKEysyYUbqxiwx1L1RpZ7XER\nibTUlnoiPqBGjVKJ4f6OLXZ8j92bkDGbgGmakUikqakpEomUeKLeQ5BcY1lWc3OzbdtBo7wm3CAi\nDZHlmXuTqcEBe9N2e6OzLo1QWaSllmOmDFVrJAQpZwrUyI4nOlesSyWGc9kETNNsbW3lmCkbQXKH\nZVn19fVLI8trw0uDRsXEg5OpwXh/15933BturGKfEqVrJAQpN2rUKHZtR16bwBb74Wg02traOj0z\n9D6uIbnAqVFL9J5KszqX8UGjoiGy/AyjorPjTtF+q1K7RsiFMjXKexPor4jFYiJCkxwEqVj51iij\nJrxURDRvEjWChjXKcDYBmpRBkIpScI0cmjeJGkHnGjloUjaCVLgia+TINMkIzdPqNQ7UCNTIkWnS\nggULNH+NA2+MLVx9fX10WVsxC9FRE166rPqbzqVvTVAjKFAjEYld2+HWJhAJR2+//XZXZuVfBKlA\nziF2pbnElUerNJekEsN2POHKo3kcNYIaNerv2CqubgLOu+ldeTSf4pRdgXp6ekRk0hd3jvL6m4Or\nHrotnU5/5dofnT7vrMztQaOi0qy2exPK79HUCGrUSEScvyBz2QQe+tPPB3Zvuu7qFvPs907LP7fr\nL9aT9w/s7jvphNmXXlTfELmx0qzu6enR+Q2zHCEVyLIs58xv7kZGRlY9+P0dLz69M7H5yJFDo+4N\nGuXOH1wKo0ZQpkYiYvcmctkE7Je2rt+wZmdi84GDb2dufHbHE//yH3//1LZHyuYEh99Orn9izcr7\nvx00yp1TL9oiSAWybXuReVlev/KHx1btSGw+3r214YZUYrjoeXkXNYJKNRKRVGJ44k1g/RO/ufu+\nW+9Y9eVRf4CmJf3rtT86evTw5z99a9vXO75343/MmDHj2R1PVC64zLbtqZ20t3HKrhAFXEDatWdL\nV8+/Xbzwo1uef3zcAUGjXETsuJpn7agRFKtRLheQ4v1de14eGHv79p1/SaYGjbln1l3xORGpmH9+\n69c6jhw9HAgERMSyLG3P2hGkQvT09FSa1blfQDp0eP+qB2+bM7vsS9fcdmv71SIigcCoMc5lpP6O\nrept2dQIitVIROx4YtJN4HNXrXBO09193y3Zt+/c84yIXHhu+JEn1jz+1O8OHzl4ycKP/bdP/N2J\nJ5xcaVavXr2aICEPzieo5j7+vnXtr76RuOn6fy479fQJhgWN8qHE00XPzluoEdSrkYikEsPlxnkT\nj1l83uXj3v7m8Csist3euHHrf5580pyDh/Y9suFee/DZW2/4ZdAo1/msHdeQCmGapvNJ8rl4atsj\nj2/6be2l14QXRyYe6XwKeLGT8xJqBCVrJCJGqCz3TWCUw0cOichb+9788mf/50+/0/ONL/xMRHYm\nNm8ZeNz5FHAX5+kvBKkQTU1NydRQjoOtJ+8XkTf3vnz3fbfe0/EPzo1rfv/j9U+sGTUymRpSaeOm\nRlC1RiISbqzKfRMYZe6c00TkzNPPufySq0Tk4oW1ofJFIrJ7aFsyNVRXV+fiPP2FU3aFcP6EGbD7\ncnmHdlrSIrJt55PZN27bsaFsTjD7lt7+LhExaxTZu6kRFK6RiBiheZLzJjDKGaedLSInzDoxc8tJ\nJ8wWkcFXXhARbS8gCUEqjPOdj/H+tbmsxYbI8sjljc7P6fTIyvu/LSLXf+qW8865JHvYdntjuLFK\njVN21Ahq10je/fbnHDeBUS5dXH/v7+8Yem3Xrj1bzjvn4leTL9ovbRWR/QffjkajOp+yI0gFikQi\n//KTe3IZuXDBey8MTadHnB8+eGHNB4LnZg8bsPtClRO95MEvqBGUr5HDrD3nmV8Wchlp3twzai9t\neHzTb+/61U0XhD5sDz579NgR8+yqV994UeQi1+fpI1xDKlBdXZ3zPcRuPWAyNaTAp31TI2hSIxEx\na0IFbwJ/s/R7kSsaZ8084dkdTxw4+PaHKq/8+ufvTKaGmpqaXJ+nj3CEVCDTNE3TzPeAPRCYcU/b\nxrG3v3MByeebODWCPjUSESM0zwiV5bIJjH3Wz5gx4/OfvvX6T9/yWjIxb+4ZJ514irMJ6HwBSThC\nKphpmq2trQN231prZZEPNWD3xTrbog80ujKxUqFG0KpGImKEyiIttcVsAgEJzA+ee9KJpzibQHd3\nt7sz9B2OkArnfJXWt1Z8V0QaIssLe5ABu689dmP0gUZf7+PUCLrVyOGcZn+sba0UvQl0d3drfngk\nBKlIRTaJGkENetbIUWSTqFE2glSsgptEjaAGnWvkKLhJ1GgUguSC7CYtMqsnvcKZTA3G+7u6rJXU\nCH5HjRzZTcprE6BG2QLpdLrUc1CEZVnNzc3O565WmtWLzMuyv7wrmRocsDe9nhrsevf6580bvuLr\nt8FSoxxZ7XERibTUlnoi7qNGo9jxROeKdanEcI6bwK5du3R+G+xYBMlNtm3btt3T02NZlmVZIuKs\ny8ybFYxQWbixyqwJ+X0Tp0a5UzVI1GhcqcRwKrHX7k3Y8T12r/Md5+/bBEzTjEajdXV1HBiNRZCm\nim3blmX19PTEYrFIS40otCVRo7woGSRqlItUYtjuTdjxRH/H1ra2NhFpbW0t9aQ8jSBNuUAg0DbY\nUupZuIYa5Uu9IFGjfLVVtLPT5oI3xiIP1AjUCFOHICFX1AjUCFOKICEn1AjUCFONIGFy1AjUCNOA\nIGES1AjUCNODIGEi1AjUCNOGIOG4qBGoEaYTQcL4qBGoEaYZQcI4qBGoEaYfQcJo1AjUCCVBkPA+\n1AjUCKVCkPAeagRqhBIiSHgHNQI1QmkRJIhQI1AjeABBAjUCNYInECTdUSNQI3gEQdIaNQI1gncQ\nJH1RI1AjeApB0hQ1AjWC1xAkHVEjUCN4EEHSDjUCNYI3ESS9UCNQI3gWQdIINQI1gpcRJF1QI1Aj\neBxB0gI1AjWC9xEk9VEjUCP4AkFSHDUCNYJfECSVUSNQI/gIQVIWNQI1gr8QJDVRI1Aj+A5BUhA1\nAjWCHxEk1VAjUCP4FEFSCjUCNYJ/ESR1UCNQI/gaQVIENQI1gt8RJBVQI1AjKIAg+R41AjWCGgiS\nv1EjUCMogyD5GDUCNYJKCJJfUSNQIyiGIPkSNQI1gnoIkv9QI1AjKIkg+Qw1AjWCqgiSn1AjUCMo\njCD5BjUCNYLaCJI/UCNQIyiPIPkANQI1gg4IktdRI1AjaIIgeRo1AjWCPgiSd1EjUCNohSB5FDUC\nNYJuCJIXUSNQI2iIIHkONQI1gp4IkrdQI1AjaIsgeQg1AjWCzgiSV1AjUCNojiB5AjUCNQIIUulR\nI1AjQAhSyVEjUCPAQZBKiRqBGgEZBKlkqBGoEZCNIJUGNQI1AkYhSCVAjUCNgLEI0nSjRqBGwLgI\n0rSiRqBGwPEQpOlDjUCNgAkQpGlCjUCNgIkRpOlAjUCNgEkRpOlAjTTX37FVqBEwmUA6nS71HBQX\nCARKPQUAJcZOm4tZpZ6AFtoGW0o9BZRG583r7N5EuLEq0lJb6rmgBJzT9anEcKkn4g+csgOminPd\nKNxYVeqJoDQyF49LPRHfIEjAlOBVDJrjpUwFIEiA+6iR5qhRYQgS4DJqpDlqVDCCBLiJGmmOGhWD\nIAGuoUaao0ZFIkiAO6iR5qhR8QgS4AJqpDlq5AqCBBSLGmmOGrmFIAFFoUaao0YuIkhA4aiR5qiR\nuwgSUCBqpDlq5DqCBBSCGmmOGk0FggTkjRppjhpNEYIE5IcaaY4aTR2CBOSBGmmOGk0pggTkihpp\njhpNNYIE5IQaaY4aTQOCBEyOGmmOGk0PggRMghppjhpNG4IETIQaaY4aTSeCBBwXNdIcNZpmBAkY\nHzXSHDWafgQJGAc10hw1KgmCBIxGjTRHjUqFIAHvQ400R41KiCAB76FGmqNGpUWQgHdQI81Ro5Ij\nSIAINdIeNfICggRQI91RI48gSNAdNdIcNfIOggStUSPNUSNPIUjQFzXSHDXyGoIETVEjzVEjDyJI\n0BE10hw18iaCBO1QI81RI88iSNALNdIcNfIyggSNUCPNUSOPI0jQBTXSHDXyPoIELVAjzVEjXyBI\nUB810hw18guCBMVRI81RIx8hSFAZNdIcNfIXggRlUSPNUSPfIUhQEzXSHDXyI4IEBVEjzVEjnyJI\nUA010hw18i+CBKVQI81RI18jSFAHNdIcNfI7ggRFUCPNUSMFECSogBppjhqpgSDB96iR5qiRMggS\n/I0aaY4aqYQgwceokeaokWIIEvyKGmmOGqmHIMGXqJHmqJGSCBL8hxppjhqpiiDBZ6iR5qiRwggS\n/IQaaY4aqY0gwTeokeaokfIIEvyBGmmOGumAIMEHqJHmqJEmCBK8jhppjhrpgyDB06iR5qiRVggS\nvIsaaY4a6YYgwaOokeaokYYIEryIGmmOGumJIMFzqJHmqJG2CBK8hRppjhrpbFapJ6As27Yty9q9\ne7eIdN68zqwNhRurSj0pr6NGmlOsRqnEsN2bSCX2ikhzc3NdXV00Gi31pDwtkE6nSz0Hpdi2vXr1\nasuyLMsSkaBRXmlWJ1NDA3afiBihMrMmFG6sUuP55i5Va2S1x0Uk0lJb6ol4nTI1SiWG+zu22PE9\ndm9CxmwCpmlGIpGmpqZIJFLiiXoPQXKNZVnNzc22bQeN8ppwg4g0RJZn7k2mBgfsTdvtjc66NEJl\nkZZajpkyVK2REKTcqFEjO57oXLEulRjOZRMwTbO1tZVjpmwEyR2WZdXX1y+NLK8NLw0aFRMPTqYG\n4/1df95xb7ixin1KlK6REKQcKFOj2LUdeW0CW+yHo9Foa2vr9MzQ+7iG5AKnRi3ReyrN6lzGB42K\nhsjyM4yKzo47RfutSu0aYVIq1SjvTaC/IhaLiQhNchCkYuVbo4ya8FIR0bxJ1EhzetYow9kEaFIG\nQSpKwTVyaN4kaqQ5zWvkoEnZCFLhiqyRI9MkIzRPq9c4UCPNUaOMTJMWLFig+WsceGNs4err66PL\n2opZiI6a8NJl1d90Ln1rghppTo0aiUjs2g63NoFIOHr77be7Miv/4gipQM4hdqW5JK/feuhPPx/Y\nvem6q1vMs993MFRpLkl1DtvxhN+fn7mgRppTpkb9HVslt01g7BN/wO57aP2/Zo85evTwi0O2ZVk6\nvz+JIBWop6dHRCZ9cWc2+6Wt6zesOXLk0IGDb4+6K2hUVJrVdq/6QaJGmlOmRiJix533vU6yCYz7\nxN+R2LwzsXnUyEqzuqenhyAhb5ZlOWd+c7H+id88v7v/mYHHjh47crwxQaO8v+NRtV/aQI00p1KN\nRMTuTUy8CUzwxB96baeIfPGa7190/hWZG9daK2OxmM4vbeAaUoFs215kXpbj4Hh/11PbHpmgRiJS\nG25IJYbdmJpHUSPNKVYjEUklhifeBCZ44g++ukNEFi5YEjQqTis7K2hUBI2K2nCDbdtTNFtf4Aip\nEPleQPrcVSuco/W777vleGOCRrmIqHoZiRppTr0a5XIB6XhP/JGRkZdft2fMmBl/6nc9Gx84euxI\n5YLqL17zfWcT0PkyEkdIhejp6ak0q3O/gLT4vMsvvaj+0ovqJxjjXEZyVrliqJHm1KuRiNjxxKSb\nwPGe+K+9uefI0cMjI8fWb1hjzJ1/7NjRrS/E22PLjbnzK83q1atXT+XEPY0jpEI4n6Dq+sMGjfKh\nxNOuP2xpUSPNKVkjEUklhsuN8wr73X0H9i4+/4qTTpz9xYbvzZ1z+s49z9yx6obX3tiz6dn1QaNc\n57N2HCEVwjRN55Pk3eV8CrjrD1tC1EhzqtZIRIxQWcGbwPnnXLLiS//ra9e3z51zuvOPCxcsEZHE\nKwPOp4C7OE9/IUiFaGpqSqaGXH/YZGpIpectNdKcwjUSkXBjVcGbwNPbH/3j4796evujmVsCEhCR\nmTNmJVNDdXV17kzRhwhSIZw/Ydw9SOrt7xIRs0aRpy410pzaNRIRIzRPCt0EXhx67sH//FnsodZ9\nB4ZFZM/LAwO7N4nIwcP7RUTbVzQIQSqM852P8f61Lj7mdntjuLFKjVN21EhzytdI3v3258I2gSsu\nufrkk+bsP/jWd+76r3f96mv/tOqGkZFjCxcsOXDwrWg0yik75C0Sibh7hDQVF6VKghppTocaOcza\ncwp72n4geO7/+MLPFlR88NDhA8/tfPLwkUMf+dCnv3rdPymzCRSMV9kVqK6urq2tbcDuy+tzFe9p\n23i8u5KpoYbGK92YWilRI83pUyMRMWtCVntvLpvA2Cf+Bed++LvLf7Vv/97hfckzjLNPOOEkEUmm\nhpqamqZqun7AEVKBTNM0TdOts3bvXEDy+XOYGmlOqxqJiBGaZ4TKitkE5pwyr/zM850aOZuAzheQ\nhCAVzDTN1tbWAbtvrbWyyIcasPtinW3RBxpdmVipUCPN6VYjETFCZZGWWhc3ge7ublcm5l+csiuc\n81Va31rxXRFpiCwv7EEG7L722I3RBxp9/TSmRprTsEYO50s1H2tbK0VvAt3d3ZofHglBKlKRTaJG\nUIC2NXIU2SRqlI0gFavgJlEjKEDzGjkKbhI1GoUguSC7SYvM6klfcpNMDcb7u7qsldQIvkaNMrKb\nlNcmQI2yBdLpdKnnoAjLspqbm53PXa00qxeZl2V/eVcyNThgb3o9Ndj17vXPmzd8xddvg6VGObLa\n4yKi3lcvUqOx7Hiic8W6VGI4x01g165dOr8NdiyC5Cbbtm3b7unpsSzLsiwRcdZlMjXkvOXNCJWF\nG6vMmpDfn8PUKHdKBokaHU8qMZxK7LV7E3Z8j93rfMf5+zYB0zSj0WhdXR0HRmMRpKli27ZlWT09\nPbFYLNJSIwptSdQoL+oFiRrlKJUYtnsTdjzR37G1ra1NRHT+evJcEKQpFwgE2gZbSj0L11CjfCkW\nJGpUgLaKdnbaXPDGWOSBGmmOGmFKESTkihppjhphqhEk5IQaaY4aYRoQJEyOGmmOGmF6ECRMghpp\njhph2hAkTIQaaY4aYToRJBwXNdIcNcI0I0gYHzXSHDXC9CNIGAc10hw1QkkQJIxGjTRHjVAqBAnv\nQ400R41QQgQJ76FGmqNGKC2ChHdQI81RI5QcQYIINdIeNYIXECRQI91RI3gEQdIdNdIcNYJ3ECSt\nUSPNUSN4CkHSFzXSHDWC1xAkTVEjzVEjeBBB0hE10hw1gjcRJO1QI81RI3gWQdILNdIcNYKXESSN\nUCPNUSN4HEHSBTXSHDWC9xEkLVAjzVEj+AJBUh810hw1gl8QJMVRI81RI/gIQVIZNdIcNYK/ECRl\nUSPNUSP4DkFSEzXSHDWCHxEkBVEjzVEj+BRBUg010hw1gn8RJKVQI81RI/gaQVIHNdIcNYLfESRF\nUCPNUSMogCCpgBppjhpBDQTJ96iR5qgRlEGQ/I0aaY4aQSUEyceokeaoERRDkPyKGmmOGkE9BMmX\nqJHmqBGURJD8hxppjhpBVQTJZ6iR5qgRFEaQ/IQaaY4aQW0EyTeokeaoEZRHkPyBGmmOGkEHBMkH\nqJHmqBE0QZC8jhppjhpBHwTJ06iR5qgRtEKQvIsaaY4aQTcEyaOokeaoETREkLyIGmmOGkFPBMlz\nqJHmqBG0RZC8hRppjhpBZ7NKPQG8R7EapRLDdm8ildhrx/cYoTKzNhRurCr1pDxNsRqxAJAvguQV\nytQolRju79hix/fYvQkRCRrllWZ18smhzo51nTevM0JlZk0o3Filxp7rImVq/PHb6wAABxZJREFU\nxAJAwQiSJ6hRIzue6FyxLpUYDhrlNeGGiyPSEFmeuTeZGhywN223NyafHIp1dBihskhLLX8yO9So\nEQsARSJIpadMjWLXdiyNLK/97NKgUTF2QNCoqAlX1ISXikgyNRjv77La700l9kZaaqd9st6iTI1Y\nACgSQSoxlWrUEr2n0qzOZXzQqGiILD/DqOjsuFNEdN6SVKoRCwBFIkilpGeNMpw/lnXekvSsUQYL\nAKMQpJLRvEYOnbckzWvk0HkBYCyCVBrUKCOzJRmhefpc4qZGGXouAIyLN8aWgBo1EpHYtR3RZW3F\nbEaOmvDSZdXftNrjrszK+9SokbAA4DaOkKabMjXq79gqIpXmkhzH9/Z3Pdr34OCrL8ydc/qli+s/\n/fEbZp98aubeSnNJqnPYjif8vkdPSpka5bUAjhw9/Lvuu59+znp7f+qC0Ic/uuQz4cWR7AH6LABM\ngCOkaaVMjUTEjjtvexznBb5j/eGxVbHOtp2JzQcP7X/tjT0Px//Pv6755sjISGZA0KioNKudt1Iq\nTJkaSZ4L4Ge//sbDj//qleSLR48d2Tzw2D0dt27c8nD2AE0WACZGkKaPSjUSEbs34Zz9n9S+A8Pr\n/hwTkYbI8p9827r+U7eIyPO7N23f9ZfsYUGj3PmjW1Uq1UjyWQD9z1kDdt+smSf8w5f//affffSa\n+q+OjIzc/8e70umR7GHKLwBMiiBNE8VqJCKpxPAi87JcRr449NzBQ/vnzC5zTtNFPtJ4Wtl8EXnp\n1Reyh9WGG1KJ4SmZqwcoViPJZwHsSGwWkQ9eWHP+OZcEJPCpK5tPOnF26q3Xnt/9VPYwtRcAckGQ\npoN6Ncrr+kFA5OKFtVd86FMzZswUkZFjxw4dOSgixtwzs4cFjXJ590SQYtSrUV4L4O39KRE5YdaJ\nzj/OmDFz1swTROSlV973F4nCCwA54kUN06G/Y2ukpUalFxHZ8T2VZnWO1w8Wn3/F4vOvcH5OS7pj\nXfv+A8OnnDx30XmXZw9zriJY7b2mWhcSnI8ZNWtCdm9CmWskeS2As84wRWTbziff2vfG3Dmnb3k+\nvu/AsIi8tf/N7GGqLgDkjiBNue7u7p6enlLPwmWxoZjz92xe3tr3Zqyzbcvzj8+aecKXPvODuXNO\nGzUgaJQPDSUjc69yaZrecJWIWv9CkucC+Oil1/z+0V/uPzD8j7+4vuLMC3Yknp45c9axY0fT6fSo\nkWouAJHW7h+Vegr+QJCmXCQSiUQipZ6Fy2zbXtv5x7x+ZWdi8z33fzs1/OppZfNv+OwPFy4Y52zP\ngN3XsOyq1tZWl6aJqZLXAjj1FONb0ZX/3tn60isv7D/w1scv++yeV54fsPtOK/vAqJEsAM0RJBSi\nqakpFovlPn77ro0//83Nh48cvHhh7d/+9Q/nzC4bd1gyNVRXV+fOFDGV8l0AofJFP7jp3rf2vTn7\npDkzZsxqueMTInL2/AtGDWMBaI4XNaAQpmmKyIDdl8vgQ4cP/O//+73DRw5eUvmxr3/+J8erUW9/\nl4iodzSppLwWwKZnH7npHz/ynbuWzj5pzqxZJz75zLr9B4bnnx664NwPZw9jAYAjJBTCNM1IJBLv\nX5vLx8Zs2PyH4beTIrLl+fjXf/hfMrf/1Uc+/7mrVmT+cbu9MRqNOjsdPC6vBXDR+Veceorxxt6X\nf/zLvz3nrMqNWx4OBGZc81c3jRrGAgBHSChQJBLJ8Q/kl1553vkhnR7JNuqado6PBo/IfQHMPvnU\nv/vvd1XMv+DFoefiT/2u7NRg02d+cPnFnxw1jAWAwNgXugC5sCyrvr6+yE96znZj22Xd3d2csfGL\nfBdAWtKp4dcknT5t3ujXMjhYAOAICQUyTdM0zXj/WlcejesHvpPvAghI4LSy+cerEQsAIjKzra2t\n1HOALxmGYRjGmvt+vf/g24uKO0gasPt+ce+3uru7uX7gIywAuI4goXDhcHj+WcFVq39RzJY0YPe1\nx27kXI0fsQDgLoKEohS5JbEZ+R0LAC4iSChWwVsSm5EaWABwC0GCC7K3pEAOX9qWTA2uf2LN6s7b\n2YzUwAKAK3jZN1xjWVZzc7Nt20GjvNKsXmRelv0FbsnU4IC96fXUYJe10rll165dXMRWCQsARSJI\ncJNt27Zt9/T0WJZlWZaIOHtTMjXkvO3RNM1oNFpXV8ffxUpiAaAYBAlTxbZty7Kcr95w/hDmU5y1\nwgJAvggSAMAT+KQGAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACAJxAk\nAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACA\nJxAkAIAnECQAgCcQJACAJxAkAIAnECQAgCcQJACAJxAkAIAnECQAgCf8f+AVijuGsxoOAAAAAElF\nTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "[node,elem] = squaremesh([0 1 0 1],0.5);\n",
    "totalEdge = sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2);\n",
    "[i,j,s] = find(sparse(totalEdge(:,2),totalEdge(:,1),1));\n",
    "edge = [j,i]; \n",
    "bdEdge = [j(s==1),i(s==1)];\n",
    "showmesh(node,elem);\n",
    "findedge(node,edge,s==1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first line collects all edges from the set of triangles and sorts the column such that `totalEdge(k,1)<totalEdge(k,2)`. The interior edges are repeated twice in `totalEdge`. We use the summation property of `sparse` command to merge the duplicated indices. The nonzero vector `s` takes value 1 (for boundary edges) or 2 (for interior edges). We then use `find` to return the nonzero indices which forms the `edge` set. We can also find the boundary edges using the subset of indices pair corresponding to the nonzero value 1. Note that we switch the order of `(i,j)` in line 3 to sort the edge set row-wise since the output of `find(sparse)` is sorted column-wise. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To construct `edge` matrix only, the above 3 line code can be further simplified to one line:\n",
    "`edge = unique(sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2),'rows');`\n",
    "The `unique` function provides more functionality which we shall explore more later. However, numerical tests show that the running time of `unique` is around 3 times of the combination `find(sparse)`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Node Star\n",
    "\n",
    "The `elem` matrix, by the definition, is a link from triangles to vertices, i.e., `elem` is `elem2node`. The link from vertices to triangles, namely given a vertex, to find all triangles containing this vertex, is stored in the sparse matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     1     0     0     1     1     0     0     0     0\n",
      "     0     1     0     0     1     1     0     0     0\n",
      "     0     0     0     1     0     0     1     1     0\n",
      "     0     0     0     0     1     0     0     1     1\n",
      "     1     1     0     0     1     0     0     0     0\n",
      "     0     1     1     0     0     1     0     0     0\n",
      "     0     0     0     1     1     0     0     1     0\n",
      "     0     0     0     0     1     1     0     0     1\n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gkHFRYi1IrfSQAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1TZXAtMjAxOCAxNDoyMjozNOymmGsAABuM\nSURBVHic7d1drBxl4cfx39KCMSivAWILnEGa3pgYiaBQoLsLJIgEKMH4T7Bx95h4RwyBEKKYdI9B\nQggqkHCjaPdEuTFYFbS+AJ7ZKlRe0lCRiIGmA4gRiFBArIGW87+YYc/pnn2Z3Z3d5+37SS9OD53p\nw+mT+e7OPLNTWlxcFAAAph1megAAAEgECQBgCYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAV\nCBIAwAoECQBgBYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAVCBIAwAoECQBgBYIEALACQQIA\nWIEgAQCsQJAAAFYgSAAAKxAkAIAVCBIAwAoECQBgBYIEALACQQIAWIEgAQCssNr0APwXx3Gr1TI9\nCgDGlMvlSqViehQOKC0uLpoeg+dKpVLl+nNMjwJm7Hvprad++kylIg5HwWo0JIkjbR4EaeJKpVLj\nn9ebHgUM2PfSW7/4vx9EkSoVbdliejQwYXZWkppNgpQL15CAiUhrtHUr743CldZo61bT43AHQQKK\nR41AjUZAkICCUSNQo9EQJKBI1AjUaGQECSgMNQI1GgdBAopBjUCNxkSQgAJQI1Cj8REkYFzUCNSo\nEAQJGAs1AjUqCkECRkeNQI0KRJCAEVEjUKNiESRgFNQI1KhwBAkYGjUCNZoEggQMhxqBGk0IQQKG\nQI1AjSaHIAF5USNQo4kiSEAu1AjUaNIIEjAYNQI1mgKCBAxAjUCNpoMgAf1QI1CjqSFIQE/UCNRo\nmggS0B01AjWaMoIEdEGNQI2mjyABnagRqJERBAk4BDUCNTKFIAFLqBGokUEECchQI1AjswgSIFEj\nUCMLECSAGoEaWYEgIXTUCNTIEgQJQaNGoEb2IEgIFzUCNbIKQUKgqBGokW0IEkJEjUCNLESQEBxq\nBGpkJ4KEsFAjUCNrESQEhBqBGtmMICEU1AjUyHIECUGgRqBG9iNI8B81AjVyAkGC56gRqJErCBJ8\nRo1AjRxCkOAtajSCJDE9gkJRI7esNj0AYCKoUX5Jovl5JYmaTUmKIlUqKpezL9xFjZxDkOAhapTT\n3JwajezrKFKjoXI561M7TvW6tmwxNsKRUSMXEST4hhrlVK0qSbIIdfys6nUliZJErZaaTSWJY0d2\nauQoggSvUKOc0hrt3dvzD0RRdsquVlO1qtlZZ47v1MhdLGqAP6hRTgNrtFwUaWFBcZwd6C1HjZxG\nkOAJapTTUDVKudIkauQ6ggQfUKOc5uYUx0s1evtt3XST1q3TkUfqjDP07W/rwIHuG7ablC52sBA1\n8gBBgvOoUX7Npur1pd9eeaVuuUWvvaZPfUp79uib39RXvtJz23TF3dzc5Ec5PGrkB4IEt1GjoSSJ\nyuXs62ef1cMP6+ij9be/6ZFHtGuXDj9cP/6xXnml5+bpovA4nspYc6NG3iBIcBg1Gkp6tq39Dumt\ntyTpxBO1Zo0kRZE+8hFJev/9nnuIIklqtSY2xOFRI58QJLiKGg2r1TrkfN2nP61TTtFzz+n22/Xc\nc5qb0xtv6Jxz9LGP9dxDuhDcnndI1MgzBAlOokYjiOOl83WSVq3Sn/+sk07SDTdo/XrdfLPWrdPv\nfjdgJ5WKLZ93R438Q5DgHmo0grQi6Tm31Msv64or9MorWrNGX/6yTj9dzz+vyy/Xvn399jMzM8lR\n5kaNvESQ4BhqNJooyj4KqO3ee/Xkkzr1VL34oubntXu3okhxrAceGLAr4++QqJGvCBJcQo3G0XH5\nJ43Tpk1atUqSjjxSGzZI0o4d/XYyP3/Ihajpo0YeI0hwBjUaU612yJub006TpEcfzZbVvfuuHntM\nkk4/vd9OOi5ETRk18htBghuo0fjS9QjtN0mbN+uII/TkkzrnHN14oz77We3Zo2OP1Re/2HMP6cJx\nU/8E1Mh7BAkOoEaFSBdtty8jnX22fvlLffKTevxx3XabnnpK55+v3/9eH/94zz2k2y5fGTE11CgE\nBAm2o0YFqtWy5xulPvc57d6tV1/V00/rjTe0Y4fOPLPntum7q/YD/aaJGgWC5yHBatSoWJWKXnhB\n1aoWFpbe6Jxwgk44YcCGSaLZWUWRgafHUqNwECTYixoVLopUq0nqbFJ/aY0kLSxMbmjdUaOgECRY\nihpNyLBNokaYGoIEG1GjiVrepPQh5St/zkmS3UjbbGYPQ5oyahQgggTrUKMpSJs0M6P5eVWriqJs\nDV65rFZLcZytDk+fgcR1I0wHQYJdqNHUpLGp17Plc61WtoIujVOjYaBDKWoULIIEi1AjI9pl2rpV\nSWLmNqM2ahQy7kOCLaiRDagRDCJIsAI1AjUCQYJ51AjUCCJIMI4agRohRZBgEjUCNUIbQYIx1AjU\nCMsRJJhBjUCN0IEgwQBqBGqElQgSpo0agRqhK4KEqaJGoEbohSBheqgRqBH6IEiYEmoEaoT+CBKm\ngRqBGmEggoSJo0agRsiDIGGyqBGoEXIiSJggagRqhPwIEiaFGoEaYSgECRNBjUCNMCyChOJRI1Aj\njIAgoWDUCNQIoyFIKBI1AjXCyAgSCkONQI0wDoKEYlAjUCOMiSChANQI1AjjI0gYFzUCNUIhCBLG\nQo1AjVAUgoTRUSNQIxSIIGFE1AjUCMUiSBgFNQI1QuEIEoZGjUCNMAkECcOxpEZJYvJvD5wlNWIO\n+Ge16QHAJWZrlCSan1eSqNmUpChSpaJyOfsC02G2RswBvxEk5GWwRnNzajSyr6NIjYbK5ezY1D4w\n1evasmXaAwuNwRoxB0JAkJCLwRpVq0qS7ADU8bfX60oSJYlaLTWbShLz55E8ZrBGzIFAlBYXF02P\nwXOlUqnxz+tNj2Isxmu0d+/gP5kkqlZVqVh3PJqbk+T8K3fjNXJ6DpRK4kibB4saMIATNZIURVpY\nUBxnh04UyIkaiTngPoKEflypUYrj0SS4UqMUc8BpBAk9mV3FEMedR6L5eW3YoKOO0rp1uuEGvflm\nlw3bx6P0QjfGZHYVw8o5IGnvXp13njZs0Isvdt+QOeAugoTuzK7wbjZVrx/ynVtuUb2unTv19tva\ns0e3367LL9fBg122TVdbpVduMA6zK7xXzgFJBw9q82Y98oh27tT+/T23ZQ44iiChC+N3vyaJyuWl\n377+um69VZIaDe3bp7vukqQdO/SHP3TfPF0QHMcTH6fHjN/92jEHUrfcokcfzbU5c8BFBAmdjNco\nPdOy/NXxrl16+20dd5xuuklHH61rrtHJJ0vS009330MUSVKrNdFh+sx4jVbOAUmPPaZvfUuf/3yu\nPTAHXESQcAjjNZLUanUeiUolXXKJvvQlrV4tSQcO6J13JGnt2u57SO/b59XxaIzXSN3mwH/+o82b\ndeyxuuee7DulUr89MAdcRJCwxIYaSYrjznM1F16o7duzM3WLi7r2Wr3xho45Rhdc0HMnlQqfdTYK\nG2qkbnPg2mv1/PP6/vd10kl5d8IccA6f1ICMJTVKjyDp+ZaVXntN9bq2b9cRR+hHP9IJJ/Tcz8zM\nBAbnO0tqtHIObNumH/5Qs7PatEnvv593P8wB5/AOCZI1NZIURdnHwKy0c6fOOEPbt+vkk/Xgg7ry\nygG74tXxUCypkbrNgbvvlqSXXtJVV+kLX8i+ec01uvPOAbtiDriFIMGiGqW6nvpfWNBFF+nll3XJ\nJdq9Wxs3DtjJ/HyXRcPoxZ4apTrmQPqxOw89pG3b9POfZ9988EE9+WS/nTAHnEOQQmdbjSTVap0v\nbN95R1dfrf/+V5deqgce0HHHDd7JyosQ6MW2GmnFHGg0dN992a+f/jT75l136Wtf67cT5oBzCFLQ\nLKyRPrgWvfwF8k9+on/9S5J+8xt96ENavTr7dX2PD61NFw1b9T9lLQtrpBVzYONGXXXV0q/UxRfr\nrLN67oE54CKCFC47a6QPFuwuv4Twl79kX7z/vg4eXPrV6/p2um2vlRFos7NG6jYHhsUccBFBCpS1\nNUrVatmzbVJ3363FxS6/vve9Ltumr6zbD3NDL9bWKNUxB9oOOyz711+/vue2zAFHEaQQWV4jSZWK\n6vXsw56HkiSanVUUOf/8oUmzvEZiDgSJ+5CCY3+NJEWRajVJqla1sJD3xEt6JJK0sDC5ofnA/hqJ\nORAkghQWJ2qUGvZ4xJEoJydqlGIOhIYgBcShGqWWH48qFdVqXUaeJNlNlM1m9iAc9OFQjVLMgaAQ\npFA4V6NUejyamdH8vKpVRVG2/qpcVqulOM5WBqfPv+GaQX/O1SjFHAhHaTG9BxoTUyqVGv/scb/M\ntDhaow7p0qn0hbC0dGCy/BiUPibO+CAdrVEHR+dAqSSOtHkQpIkzHiQ/atQhSZy5xcSGIPlRow4O\nzQGClBPLvj3nZY3EDY/D8LJGYg74iCD5zNcaIT9fawQvESRvUSNQI7iFIPmJGoEawTkEyUPUCNQI\nLiJIvqFGoEZwFEHyCjUCNYK7CJI/qBGoEZxGkDxBjUCN4DqC5ANqBGoEDxAk51EjUCP4gSC5jRqB\nGsEbBMlh1AjUCD4hSK6iRqBG8AxBchI1AjWCfwiSe6gRqBG8RJAcQ41AjeArguQSagRqBI8RJGdQ\nI1Aj+I0guYEagRrBewTJAdQI1AghIEi2o0agRggEQbIaNQI1QjgIkr2oEagRgkKQLEWNQI0QGoJk\nI2oEaoQAESTrUCNQI4SJINmFGoEaIVgEySLUCNQIISNItqBGoEYIHEGyAjUCNQIIknnUCNQIEEEy\njhqBGgEpgmQSNQI1AtoIkjHUCNQIWI4gmUGNQI2ADgTJAGoEagSsRJCmjRqBGgFdEaSpokagRkAv\nBGl6qBGoEdAHQZoSagRqBPRHkKaBGoEaAQMRpGmgRoFrNiVqBAxSWlxcND0Gz5VKJdNDAGAYR9o8\nVpseQBCYisGanVUcq17Xli2mhwITkkRnbDxq30tvmR6IGzhlB0xKet2oXjc8DJiS1mjT9z5neiDO\nIEjARLCKIXDtGkUbTjE9FmcQJKB41Chw1Gg0BAkoGDUKHDUaGUECikSNAkeNxkGQgMJQo8BRozER\nJKAY1Chw1Gh8BAkoADUKHDUqBEECxkWNAkeNikKQgLFQo8BRowIRJGB01Chw1KhYBAkYETUKHDUq\nHEECRkGNAkeNJoEgAUOjRoGjRhNCkIDhUKPAUaPJIUjAEKhR4KjRRBEkIC9qFDhqNGkECciFGgWO\nGk0BQQIGo0aBo0bTQZCAAahR4KjR1BAkoB9qFDhqNE0ECeiJGgWOGk0ZQQK6o0aBo0bTR5CALqhR\n4KiREQQJ6ESNAkeNTCFIwCGoUeCokUEECVhCjQJHjcwiSECGGgWOGhlHkACJGgWPGtmAIAHUKHTU\nyBIECaGjRoGjRvYgSAgaNQocNbIKQUK4qFHgqJFtCBICRY0CR40sRJAQImoUOGpkJ4KE4FCjwFEj\naxEkhIUaBY4a2YwgISDUKHDUyHIECaGgRoGjRvYjSAgCNQocNXICQYL/qFHgqJErCBI8R40CR40c\nQpCwJElMj6Bo1GhYns0BauSW1aYHAMOSRPPzShI1m5IURapUVC5nXziNGuXk6xygRs4hSOGam1Oj\nkX0dRWo0VC5nx6b2gale15YtxkY4DmqUh8dzgBq5iCAFqlpVkmQHoI5XwfW6kkRJolZLzaaSxL3D\nOjXKw+M5QI0cVVpcXDQ9Bs+VSiXbfsbpkWjv3sF/MklUrapScel4ZFuN5uYkWfcmw+M5YGGNGmu+\nw5E2DxY1BCf/kUhSFGlhQXGcHeXtZ1uN7OTxHLCwRsiPIIVlqCNRyqHjETXKw+M5QI1cR5ACMjen\nOB7uSJRqH4/SC912okZ5eDwHqJEHWNQQkGZT9foh34ljfeMbnX/s0Ue7bJuutpqb69yDJahRTsvn\nwPbtuvnmLn/mjjv0mc90+b7Nc4Aa+YEgBSRJVC4f8p2dO7VzZ97Ny2U1Gopj6+5NoUb5LZ8Dr77a\n/V//zTd7bm7nHKBG3iBIoUjPtHS8tn3mGUn6wQ900UWD9xBFktRq2XUwokb5dcyB88475Of297/r\n1lt1/PH6xCd67sHCOUCNfEKQQtFqdTnTkgZp40ZFkQ4e1KpV/faQ3rcfxxatYKZGQ+mYA+vWad26\npd9edpkOO0z33qs1a3ruwbY5QI08w6KGUMRx5/m6gwf17LNavVpbt+qYY/TRj+qSS/SPf/TbSaVi\n0WedUaNhrZwDbQ8+qF/9Sl/9qi6+eMBO7JkD1Mg/BCkI6REkPd/StmeP/vc/HTigO+/U2rV67z39\n9reqVPTeez33MzMzyVEOgxoNq+scSB08qOuuU6mk664bvB9L5gA18hJBCkIUZR8Ds9zrr+vCC3XF\nFUoSPfOM/vhHrVqlPXt03339dmXDq2NqNIKucyB1zz3661916aVavz7XrozPAWrkK4IUivTU/3Jn\nn62HHtIvfqETT8x+u3GjJO3e3XMn8/Pml/xSo5GtnAOp9Ie5aVOunRifA9TIYwQpFLVa5wvb++/X\nbbfp/vuXvlMqSdLhh/fcSZ+LENNBjcaxcg5I2r9fu3ZJ0gUX5NqJ2TlAjfxGkEKRXote/gJ51y7d\neKNqNb3+uiTt3q0dOyTprLO67yFdNGxwvS81GtPKOSDpiSf03ntau1annTZ4D2bnADXyHkEKRbpg\nd/klhKuv1lFHad8+nXqqLrpI556rAwe0caMuu6z7HtJtu14VnwJqNL6Vc0DSI49k/ykPg3OAGoWA\nIAWkVsuebZNav17bt+vMM/XOO3r4Ye3fr82b9bOfZSfuOqSvrNsPc5syalSUjjkg6U9/kqS1awdv\na3AOUKNAcGNsQCoVvfCCqlUtLGQvcs89V088oX//W6+8otNO04c/3H3DJNHsrKLIzO2Q1KhAK+fA\nr3+da0ODc4AahYMgBSSKVKtJOuR4JOn443X88T23So9EkhYWJj7ClahRsXrNgf4MzgFqFBSCFJZh\nj0fUyD8OzQFqFBqCFJzlx6NKRbVal0VTSZLdRNlsZg/CmT5qNDlOzAFqFCCCFKL0eDQzo/l5VauK\nomz9VbmsVktxnK0MTp9/w3UjL1k+B6hRmEqLi4umx+C5Uqlk8884XTqVvhCWlg5MBj/O2bMazc1J\nsuXjsbuybQ74V6PGmu9wpM2DIE2c5UFaLkmM3WbU5lmN5EKQljM+B/yrkQhSbtyHhCXUCNQIBhEk\n2IIaBY4agSDBCtQocNQIIkiwATUKHDVCiiDBMGoUOGqENoIEk6hR4KgRliNIMIYaBY4aoQNBghnU\nKHDUCCsRJBhAjQJHjdAVQcK0UaPAUSP0QpAwVdQocNQIfRAkTA81Chw1Qn8ECVNCjQJHjTAQQcI0\nUKPAUSPkQZAwcdQocNQIOREkTBY1Chw1Qn4ECRNEjQJHjTAUgoRJoUaBo0YYFkHCRFCjwFEjjIAg\noXjUKHDUCKMhSCgYNQocNcLICBKKRI0CR40wDoKEwlCjwFEjjIkgoRjUKHDUCOMjSCgANQocNUIh\nCBLGRY0CR41QFIKEsVCjwFEjFIggYXTUKHDUCMUiSBgRNQocNULhCBJGQY0CR40wCQQJQ7OhRkli\n8m8PHDXChKw2PQA4xmCNkkTz80oSNZuSFEWqVFQuZ19gOqgRJocgYQimajQ3p0Yj+zqK1GioXM76\n1I5Tva4tW6Y9sNBQI0wUQUJepmpUrSpJsgh1vBOq15UkShK1Wmo2lSRc1pogaoRJI0jIxWyN9u7t\n+QeiKDtlV6upWtXsLE2aCGqEKWBRAwaztkbLRZEWFhTH2WhRIGqE6SBIGMCJGqVo0iRQI0wNQUI/\nBlcxxPFwNUq1m5QudsCYqBGmiSChJ4MrvJtN1evd/9PXv64NG/T44z23TVfczc1NZmQhoUaYMoKE\n7sze/ZokKpe7fP/xx3Xnndq5U2++2W/zdFF4HE9mcGGgRpg+goQuzNYoPdvW8Q7pjjt01VU6/3zt\n3z94D1EkSa1W0SMLBjWCEQQJnYx/MlCr1eV8XbOpbdv07ru59pAuBOcd0mioEUwhSDiE8RpJiuMu\n5+u++11t26Zt2/LupFLh8+5GQY1gEDfGYokNNUorkp5zW+6CC4bbz8xMEaMJDDWCWbxDQsaGGkmK\nouyjgMbHO6ShUCMYR5AgWVOjVCGXf+bney4cx0rUCDYgSLCrRpJqtQLe3HS9EIWuqBEsQZBCZ1uN\n9MF6hHHeJKULx3lIUh7UCPYgSEGzsEb6YNH2OJeR0m1XroxAB2oEqxCkcNlZo1Stlj3faATpu6v2\nA/3QCzWCbVj2HSibaySpUtELL6ha1cJC5xudxcV+GyaJZmcVRTw9dgBqBAsRpBBZXiNJUaRaTVL3\nJvWS1kjSwsLkhuYDagQ7EaTg2F+j1LBNokY5USNYiyCFxZUapZY3KX1I+cqFc0mS3UjbbGYPQ0If\n1Ag2I0gBcatGqbRJMzOan1e1qijK1uCVy2q1FMfZ6vD0GUhcN+qPGsFyBCkULtYolcamXs+Wz7Va\n2Qq6NE6NBh3KhRrBfgQpCO7WaLl2mbZuVZJwm9EQqBGcwH1I/vOjRh2oUX7UCK4gSJ7zskbIjxrB\nIQTJZ9QocNQIbiFI3qJGgaNGcA5B8hM1Chw1gosIkoeoUeCoERxFkHxDjQJHjeAuguQVahQ4agSn\nESR/UKPAUSO4jiB5ghoFjhrBAwTJB9QocNQIfiBIzqNGgaNG8AZBchs1Chw1gk8IksOoUeCoETxD\nkFxFjQJHjeAfguQkahQ4agQvEST3UKPAUSP4iiA5hhoFjhrBYwTJJdQocNQIfiNIzqBGgaNG8B5B\ncgM1Chw1QggIkgOoUeCoEQJBkGxHjQJHjRAOgmQ1ahQ4aoSgECR7UaPAUSOEhiBZihoFjhohQATJ\nRtQocNQIYSJI1qFGgaNGCBZBsgs1Chw1QsgIkkWoUeCoEQJHkGxBjQJHjQCCZAVqFDhqBIgg2YAa\nBY4aASmCZBg1Chw1AtoIkknUKHDUCFiOIBlDjQJHjYAOBMkMahQ4agSsRJAMoEaBo0ZAVwRp2qhR\n4KgR0AtBmipqFDhqBPRBkKaHGgWOGgH9EaQpoUaBo0bAQARpGqhR4KgRkEdpcXHR9Bg8VyqVJDUa\npscBQ+JYcazonFOiDSebHgvMiL+zkyNtHgRp4uI4brVapkcBwJhyuVypVEyPwgEECQBgBa4hAQCs\nQJAAAFYgSAAAKxAkAIAVCBIAwAoECQBgBYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAVCBIA\nwAoECQBgBYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAVCBIAwAoECQBgBYIEALACQQIAWIEg\nAQCsQJAAAFYgSAAAKxAkAIAVCBIAwAoECQBgBYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAV\nCBIAwAoECQBgBYIEALACQQIAWIEgAQCsQJAAAFYgSAAAKxAkAIAV/h9l2LLo/iwqggAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "NT = size(elem,1); N = size(node,1);\n",
    "t2v = sparse([1:NT,1:NT,1:NT], elem, 1, NT, N);\n",
    "display(full(t2v));\n",
    "nodeStar = find(t2v(:,5));\n",
    "showmesh(node,elem);\n",
    "findnode(node,5);\n",
    "findelem(node,elem,nodeStar);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NT x N matrix `t2v` is the incidence matrix between triangles and vertices. `t2v(t,i)=1` means the i-th node is a vertex of triangle t. If we look at `t2v` column-wise, the nonzero in the i-th column of `t2v(:,i)` will give all triangles containing the $i$-th node. Since sparse matrix is stored column-wise, the star of the $i$-th node can be efficiently found by `nodeStar = find(t2v(:,i))`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We label three edges of a triangle such that the i-th edge is opposite to the i-th vertex. We define the matrix `elem2edge` as the map of local index of edges in each triangle to its global index. The following 3 line code will construct `elem2edge` using more output from `unique` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "totalEdge = sort([elem(:,[2,3]); elem(:,[3,1]); elem(:,[1,2])],2);\n",
    "[edge, i2, j] = unique(totalEdge,'rows','legacy');\n",
    "elem2edge = reshape(j,NT,3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Line 1 collects all edges element-wise. The size of `totalEdge` is thus 3NT x 2. By the construction, there is a natural index mapping from `totalEdge` to `elem`. In line 2, we apply `unique` function to obtain the edge matrix. The output index vectors `i2` and `j` contain the index mapping between `edge` and `totalEdge`. Here `i2` is a NE x 1 vector to index the last (2-nd in our case) occurrence of each unique value in `totalEdge` such that `edge = totalEdge(i2,:)`, while `j` is a 3NT x 1 vector such that `totalEdge = edge(j,:)`. (Try `help unique` in MATLAB to learn more examples. `legacy` is used since the version change of MATLAB.) Then using the natural index mapping from `totalEdge` to `elem`, we reshape the 3NT x 1 vector `j` to a NT x 3 matrix which is `elem2edge`.\n",
    "\n",
    "We then define a NE x 4 matrix `edge2elem` such that `edge2elem(k,1)` and `edge2elem(k,2)` are two triangles sharing the k-th edge for an interior edge. If the k-th edge is on the boundary, then we set `edge2elem(k,1) = edge2elem(k,2)`. Furthermore, we shall record the local indices in `edge2elem(k,3:4)` such that `elem2edge(edge2elem(k,1),edge2elem(k,3))=k`. Similarly `edge2elem(k,4)` is the local index of k-th edge in `edge2elem(k,2)`. \n",
    "\n",
    "To construct `edge2elem` matrix, we need to find out the index map from `edge` to `elem`. The following code is a continuation of the code constructing `elem2edge`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "i1(j(3*NT:-1:1)) = 3*NT:-1:1; i1=i1';\n",
    "k1 = ceil(i1/NT); t1 = i1 - NT*(k1-1);\n",
    "k2 = ceil(i2/NT); t2 = i2 - NT*(k2-1);\n",
    "edge2elem = [t1,t2,k1,k2];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code in line 1 uses `j` to find the first occurrence of each unique edge in the `totalEdge`. In MATLAB, when assign values using an index vector with duplication, the value at the repeated index will be the last one assigned to this location. Obvious `j` contains duplication of edge indices. For example, `j(1)=j(2)=4` which means `totalEdge(1,:)=totalEdge(2,:)=edge(4,:)`. We reverse the order of `j` such that `i1(4)=1` which is the first occurrence.\n",
    "\n",
    "Using the natural index mapping from `totalEdge` to `elem`, for an index `i` between `1:N`, the formula `k=ceil(i/NT)` computes the local index of i-th edge, and `t=i-NT*(k-1)` is the global index of the triangle which `totalEdge(i,:)` belongs to. The `edge2elem` is just composed by `t1,t2,k1` and `k2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "30px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
